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Preface 


This  study  investigated  the  state  observability 
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equations . 
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ABSTRACT 

The  relative  position  determination  of  a  cluster  of 
satellites  operating  in  a  low  earth  orbit  is  investigated. 

A  U-D  Covariance  Factorization  Kalman  Filter  is  used  for  the 
on-board  estimator  with  dynamics  based  on  the 
Clohessy-Wiltshire  equations.  Measurements  consist  of  range 
data  between  a  single  host  satellite  and  the  remaining 
cluster.  Therefore  only  relative  position  and  velocity 
states  with  respect  to  the  host  satellite  can  be  determined. 
A  15-sample  Monte  Carlo  simulation  was  conducted  with 
clusters  of  2,  5  and  10  satellites,  respectively. 

Performance  results  consist  of  average  error,  average  true 
error  and  filter  covariance  as  a  function  of  time. 
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AUTONOMOUS  NAVIGATION 


OF  A 

SATELLITE  CLUSTER 


I.  Introduction 


The  concept  of  using  a  recursive  filter  for  relative 
position  determination  for  a  cluster  of  satellites  acting 
as  a  space  based  radar  was  the  topic  of  two  previous  Air 
Force  Institute  of  Technology  masters'  theses.  The  initial 
concept  called  for  a  cluster  of  up  to  ten  satellites 
orbiting  in  a  near  circular,  low  earth  orbit.  The  cluster 
was  placed  within  a  volume  of  space  of  dimensions  500  X  500 
X  500  meters3  and  the  accuracy  requirement  of  the  filter  was 
25  meters.  The  accuracy  was  based  upon  the  requirement  to 
form  a  clear,  cohesive  image  and  is  a  function  of  the 
radar's  wavelength.  The  filter  operated  on  range  data 
determined  from  synchronized  clock  pulses. 

Captain  Michael  L.  P.  Ward  investigated  the  feasibility 
of  using  a  recursive  filter  to  determine  the  relative 
position  of  each  satellite  within  the  cluster.  The  filter 
was  a  U-D  covariance  factorization  Kalman  filter.  During 
testing.  Captain  Ward  discovered  that  the  downrange 
component  of  the  state  was  unobservable  (6:2-17).  The 
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cluster's  state  was  therefore  modified  to  include  only 
relative  downrange  components.  The  filter's  performance 
satisfied  the  required  accuracy  and  proved  promising  under 
initial  testing  (6:4-1).  The  continued  testing  of  the 
filter  was  conducted  by  Captain  Sherrie  Norton  Filer. 

Captain  Filer  investigated  the  discrepancy  between  the 
results  from  each  satellite  (3:2).  Captain  Ward's  initial 
testing  revealed  that  filter  performance  was  different  for 
each  satellite.  The  satellite  that  contains  the  filter 
under  investigation  is  defined  as  the  host  satellite. 

During  continued  research.  Captain  Filer  discovered  that 
the  filter  for  the  host  satellite,  designated  satellite  #1, 
was  unable  to  update  its  state  vector  components . 

Additional  information  was  sent  to  satellite  #l's  filter  in 
an  effort  to  improve  its  ability  to  update  its  own  state. 
When  this  effort  proved  unsuccessful.  Captain  Filer 
initiated  a  search  for  other  unobservable  state  components . 
The  research  yielded  mixed  results;  it  seems  that  none  of 
the  states  were  truly  unobservable  (3:52). 

The  purpose  of  this  thesis  is  to  continue  the 
investigation  into  the  unobservability  of  the  state 
components.  The  orbital  altitude,  cluster  radius,  and 
accuracy  requirement  will  remain  the  same.  The  cluster 
geometry  will  be  a  random  distribution  of  satellites  within 
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the  cluster  volume.  Satellite  #1  will  be  considered  the 
primary  satellite  and  its  filter's  performance  will  be  the 
source  of  the  results  presented  in  this  thesis.  Extensive 
usage  of  Captains  Ward  and  Filer's  computer  code  will  be 
made  in  an  effort  to  reduce  development  and  testing  time. 
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II.  Background 


The  U-D  covariance  factorization  filter  computer  code 
was  originally  generated  by  Captain  Ward.  Additionally,  a 
truth  model  was  developed  to  provide  the  estimator  with 
perfect  or  corrupted  range  data  and  to  provide  a  true  state 
for  comparison  against  the  estimator's  state.  The 
following  sections  detail  the  development  of  the  truth 
model  and  the  estimator  as  previously  presented  by  Captain 
Ward  (6:2-1-2-21) . 

2 . 1  Truth  Model 


The  truth  model  was  base ‘  on  two-body  orbital  dynamics 
and  provided  the  cluster's  true  state  and  the  relative 
distance  between  each  satellite.  The  range  data  is  a 
relative  measurement  and  therefore  absolute  positions 
cannot  be  determined  by  the  estimation  filter. 
Alternatively,  the  relative  distance  is  defined  with 
respect  to  a  rotating  reference  point  located  in  a  circular 
orbit  of  radius  R  with  velocity  -v<  ^  /  /? .  The  rotating 
reference  frame  (.v,y,£)  rotates  with  respect  to  an 
earth-centered  inertial  reference  frame  ( i,],k )  at  a 
constant  angular  rate  of  uok  (see  Figure  1). 
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The  position  of  each  satellite  within  the  cluster  was 
defined  randomly  about  the  reference  point  and  takes  the 
form: 


fR\ 

rt  =  <  0  >  +  {n^}*  (500meters)  (3) 

where  n;  is  a  vector  of  random  numbers  between  -0.5  and  0.5 

from  a  uniformly  distributed  random  number  generator. 

The  radial  and  out-of-plane  components  of  a  satellite's 
velocity  are  determined  by  the  Clohessy-Wiltshire  equations 
(these  equations  will  be  introduced  in  the  next  section) 
and  take  the  form:  (8:80) 

vJ=T](rJ-rre,i)  (4) 

u,£  =  r|r,£  (5) 


where  r)  is  the  mean  motion  of  the  reference  point. 


The  third  component  of  the  satellite's  velocity  is 
determined  from  the  constraint  that  the  orbital  periods  of 
all  satellites  within  the  cluster  must  be  equal  to  ensure 
that  the  cluster  remains  intact.  From  Kepler's  laws,  the 
orbital  period  of  an  elliptical  orbit  is  a  function  of  the 
semi-major  axis  of  the  orbit  (2:33).  Therefore,  the 
semi-major  axis  of  each  satellite  must  be  the  same  and 
equal  to  the  radius  of  the  reference  point.  Utilizing  the 
energy  equation,  one  may  solve  for  the  third  velocity 
component  to  obtain: 


vj 


-vj-i >li-vik'Vlk  ] 


(6) 


The  above  information  may  be  used  to  form  the  inputs  to 
the  estimator.  Captain  Ward  defined  the  initial  true  state 
as  the  position  and  velocity  components  of  each  satellite 
expressed  with  respect  to  the  rotating  reference  frame 
[ROT].  Initially,  the  axes  of  both  the  fixed,  inertial 
frame  [FIX]  and  the  rotating  frame  are  aligned  and  the 
position  and  velocity  may  be  expressed  as: 


‘[ROT] 


“(r.-fr./). 


FIX1 


(7) 
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v 

1  [ROT] 


(vv 


<■«/  )[FIX]  ' 


uo  X  r 


‘[ROT] 


(8) 


therefore  the  state  at  t=0  is 


r 1 [ROT] 
V 1 [ROT] 


X,(0)  = 


(9) 


r*[ROT) 
_  V  s[ROT) 


Additionally,  the  truth  model  outputs  an  s-1  relative 
measurement  vector. 


2,= 


I  r,-r: 


ri  ~rs  I 


U, 


(10) 


where  u,  represents  zero-mean,  white  Gaussian  noise  with  an 

associated  covariance  of  R,  (4:330).  The  noise  u,  is  the 
best  representation  for  errors  in  computing  the  range 
measurements  from  the  clock  pulses.  There  are  numerous 
sources  for  the  errors,  but  for  the  purpose  of  this  thesis 
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the  errors  are  lumped  together  into  a  single  term  and 
considered  independent  from  measurement  to  measurement 
(6:2-6) . 

The  future  position  and  velocity  vectors  for  each 
satellite  are  determined  by  the  solution  of  the  Kepler 
problem  and  the  /  and  g  equations  defined  in  terms  of  the 
eccentric  anomaly  E  (2:219).  The  equations  take  the 
following  form: 


/  =  1  —  —  (  1  —  cos  AE) 

a 


g  =  t-  A  /  — (AE-sinAE) 


/- 


\f[ias\uAE 


rr. 


a 


g  =  1  —  —  ( 1  —  cos  AE) 


where  AE  and  r  are  defined  as: 


AE  =  E/-E0 
r  =  a(  1  -  ecosEy) 


(11) 

(12) 

(13) 

(14) 


(15) 

(16) 


The  value  of  the  eccentric  anomaly,  E,,  for  any  time  t 


was  determined  by  a  Newton  iteration  scheme  due  to  the 
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transcendental  nature  of  the  Kepler  equation.  Once 
/,g,/,andg  are  determined,  the  new  position  and  velocity- 
vectors  for  each  satellite  may  be  calculated  using: 


r(0  =  f  r  0  +  gv0 

(17) 

K0  =  frB  +  gv0 

(18) 

Once  the  position  and  velocity  of  each  satellite  with 
respect  to  the  inertial  reference  frame  have  been 
determined,  the  relative  position  and  velocity  vectors  with 
respect  to  the  rotating  reference  frame  can  be  determined. 
First,  the  new  inertial  position  and  velocity  vectors  of 
the  reference  point  must  be  determined.  This  is  easily 
accomplished  by  a  rotation  about  the  k  axis  through  an 
angle  0  defined  as: 

e  =  u >-t  (19) 

Therefore,  the  new  inertial  position  and  velocity  vectors 
for  the  reference  point  are: 

R  cos  0 

rr*fUix](0=  R  sin  0  (20) 

L  o 

sin  0 

:os0 
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The  inertial,  relative  position  and  velocity  vectors 
between  each  satellite  and  the  reference  point  are 
determined  by  subtracting  the  two  solutions. 

rrel(/,x]=[r(0-rra/(0]  =  r1i>r2;  +  r3fc  (22) 

Since  the  inertial  and  the  rotating  reference  frames 

are  no  longer  aligned,  eqns  (22)  and  (23)  must  undergo  a 

coordinate  transformation  through  the  angle  0  as  previously 

defined.  The  relative  position  and  velocity  vectors 

expressed  with  respect  to  the  rotating  reference  frame  are: 

r ,  cos0  +  r2sin  G 
-r,sin0  +  r2cos0 

ri[ROT]  _  r  3 

^ / [hot j  J  v>  cos0  +  u2sin  6  +  u)  •  r2 

-vx  sin  6  +  d2cos0-  uo  ■  r , 

^3 

Once  the  position  vectors  expressed  in  the  rotating 
reference  frame  have  been  determined,  the  filter  inputs  can 
be  calculated  as  before. 

2.2  Kalman  Filter 

The  extended  Kalman  filter  is  used  due  to  the  nonlinear 
nature  of  the  measurements.  The  'U-D  covariance 
factorization'  version  of  the  Kalman  filter  was  used  to 
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solve  numerical  problems  encountered  by  Captain  Ward 
(6:2-14).  This  form  of  the  Kalman  filter  will  continue  to 
be  used  for  the  estimator  once  the  observability  problems 
are  identified  and  corrected. 

The  system's  state  propagation  is  described  in  terms  of 
a  linear  stochastic  differential  equation.  The  available 
filter  inputs  are  discrete-time,  noise  corrupted,  nonlinear 
measurements  of  the  range  between  the  satellites.  The 
Kalman  propagation  equations  for  the  system  are  (4:220): 

(25) 

Where 

Qd  =  The  covariance  of  the  dynamics  driving  noise. 

x(0=  The  estimated  state  after  time  propagation. 
x(C)=  The  estimated  state  after  the  measurement  update. 
Gd  =  Equals  the  identity  matrix  (because  the  model  is 

an  equivalent  discrete-time  representation  of  a 
continuous-time  system  (4:377)) 

4>  =  The  state  transition  matrix. 

The  time  argument  ( t‘ )  will  be  replaced  by  (±)  for  the 

remainder  of  the  text.  The  first  time  propagation  occurs 
before  any  measurement  updates,  thus  the  initial  filter 
state  and  covariance  must  be  established  (6:2-10).  For 
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ease  in  studying  the  steady  state  behavior  of  the  filter, 
x0  will  be  set  equal  to  x,(0).  The  initial  state  covariance 
matrix  P0  will  be  diagonal  with  position  elements  of  order 
10‘6A:m2  and  velocity  elements  of  order  10*12 km2 / sec2 
respectively. 

The  extended  Kalman  filter  update  equations  are  (5:44) 

*(+)  =  *(-)K<z-h[x(-)]>  (27) 

P(  +  )  =  P(-)-KHP(-)  (28) 

K  is  the  Kalman  filter  gain;  an  expression  defining  the 
gain  will  be  developed  shortly. 

The  measurement  vector  h(x,0  is  the  filter's  estimate 

of  the  range  between  the  satellites.  Therefore,  the 
residual  is  the  difference  between  the  observed  data  vector 
and  the  measurement  vector,  z-h[X(-)].  The  measurement 
vector's  form  is: 


The  matrix  H  is  developed  from  linearizing  the  h  vector 
with  respect  to  the  state  components  and  evaluating  it 
after  the  time  propagation  of  the  state,  x(-): 


13 


x-S(-) 


(30) 


H=r 


Which  has  the  following  form: 


H  = 


H,  -H,  0 

H,  0  -  H- 


0 

0 


Hs- 1  0  0  ...  -  Hs_, 


where 


Ht  = 


*i  -*2  y i ~ y 2  -z2 

hi  hi  hi 


0  0  0 


H 


s-  1 


*i-*»  yi-y*  z\-za 


/is_i  h  s- 1  h  s  ~  i 


0  0  0 


2.2.1  Dynamics 


(31) 


(32) 

(33) 


The  system  dynamics,  i.e.  the  state  transition  matrix, 
is  based  upon  the  Clohessy-Wiltshire  equations  of  motion. 
These  equations  describe  "the  relative  motion  of  two 
satellites  when  one  is  in  a  circular  orbit."  (8:78): 


x-2riy-3r|2x  =  0 

(34) 

y  +  2r\x  =  0 

(35) 

z  +  r\2  z  =  0 

(36) 
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One  may  integrate  these  equations  about  the  initial 
conditions  x0 ,  xa ,  y  0 ,  yD ,  z0 ,  and  z0  to  obtain  the  position  and 
velocity  solutions  (8:79-81) 


(2  \  x  0  2 

x(t)  =  -  -yo+3x0  cost]!  +  —  sinqi  +  4x0  +  -y0  (37) 

v  'n  /  r)  ti 


( 4y0  \  2x0  2x  0 

y(0  =  y0  -  (3y0  +  6n *.)*  +  +  6x0Jsm +  — ^costij  -  -—?■ 

(38) 


z(0  =  2  cos Tji  +  — sinrjf 

T) 


(39) 


X(t )  =  (2y0  +  3q.v0)sin  T)^  +  x0cosr]f 


(40) 


y(0  =  ~3y o  -  6tix0  +  (6t).v0+  4y0)cosr|l  -  2x0sinT^  (41) 


2(1)  =  -20risinril  +  20cost|f 


(42) 


Now,  one  may  develop  the  state  transition  matrix  (8:81). 


4-3cosip  0  0 


sin  ip 
T| 


( 1  -  cosiy)  0 


2  4  3i^ 

6(sin -  ty)  1  0  -(cosip-l)  -simp - 0 


<t>  = 


0  cosip 


simp 


3r|simp  0  0 

6r](cosip  -  1)0  0 

0  0  -  r|Sin  ip 


cosip  2simp  0 

-2simp  -3  +  4cosip  0 
0  0  cosip 


(43) 
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where  is  r)6 1  and  6 1  is  the  sample  time.  Generalizing  the 
state  propagation  equation  for  the  s  satellites  yields 


*(-) 


0  0 

0  4>2  0  •  • 

0  0  4>3  •  • 


0 


0 


0 


(44) 


2.2.2  U-D  Covariance  Factorization  Filter 


The  U-D  filter  is  especially  valuable  for  small  word 
length  micro-processors,  since  it  achieves  twice  the 
numerical  precision  capability  for  the  same  wordlength 
(4:400).  The  basis  of  the  filter  is  the  factorization  of 
the  covariance  matrices  into  a  unitary  upper  triangular  and 
a  diagonal  matrix,  such  that 


P(-)  =  U(-)D(-)Ut(-) 

(45) 

P(  +  )  =  LI(  +  )D(  +  )Ut(  +  ) 

(46) 

The  algorithm  is  initiated  with  the  same  initial 
covariance  values  as  stated  above.  First,  the  initial 
n-by-n  covariance  matrix  P0  (where  n  is  the  number  of 
states)  is  factored  into  the  UDl!T  form  by  the  following 
steps : 
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First,  for  the  nth  column  determine 


£>„„  =  P, 


(47) 


Uin  =  { 


i  =  n 

i  =  n  -  l  ,n- 2 . 1 


(48) 


Then  for  the  remaining  columns,  j  =  n- l ,  n- 2,  . . . ,  l  determine 


D„  ~  pj,- 1  D«u% 


(49) 


fl 

Pi-  I  DkkUlkUJk  /Dn  i-j-l.j-2 . 1 


(50) 


Once  the  initial  U  and  D  matrices  are  determined,  the 
state  can  be  propagated  forward  to  the  first  update  time. 
The  rc-by-2n  matrix  Y(-)  is  formed  by  augmenting  the  state 
propagation  equation  and  the  identity  matrix  G d. 


Y(-)  =  [<f>U(  +  )  |  Gd] 


(51) 


Finally,  the  2n-by-2n  matrix  D(-)  is  formed  by  using  the 
D(+)andQd  matrices  as  the  block  diagonal  elements. 


D(-)  = 


D  ( +  )  0 

o  Qd 


(52) 


The  transpose  of  the  Y(-)  matrix  forms  a  matrix  of  column 
vectors  a,  of  length  2 n. 


Y  (-)  =  [  a,  a o 


(53) 
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The  propagation  is  accomplished  by  calculating  the 
following  relationships  for  k  =  n,n~  1 . 1: 


ck  =  D(-)aJk  (ct,  =  Dy;(-)a;t,  y  =  1,2 . 2  n) 


D  u(-) 

=  aTct 

d* 

=  Ct/Dtik(-) 

(54) 

•  aJd* 

j-  1.2, 

...,k-  1 

a, 

<-  a rujk(.~)^k 

7=1.2. 

....  /c  -  1 

The  final  step  of  the  U-D  covariance  factorization 
algorithm  is  the  scalar  measurement  update.  The  following 
equations  complete  the  update  using  the  previously  computed 
UT(-)  and  D(-)  matrices,  the  1-by-n  rows  of  the  H (1,)  matrix 
and  the  measurement  covariance  value. 

f  =  UtHt 

=  D „(-)/,  j  -  1,2 . n  (55) 

a0  =  R 

Then,  for  k  =  1,  2,  ...,  n 

a  k  ~  a  k-  l  +  f  kv  k 
Dkk(  +  )  =  Dkk(-)cik_i/ak 

bk  <-  vk  (56) 

P  k  ~  ~  f  k  /  a  k-  1 

Lrjk(  +  )  =  l'lk(-)  +  bJPk  y=1.2 k-  1 

b ,  b  I  +  u  7*  1  -2 A:  -  1 

The  filter  gain  is  then  determined  by 


18 


Finally,  the  state  vector,  x(+)  and  the  covariance  matrix, 
P(+)  can  be  determined  using  Eqns  27  and  28. 

2.3  State  Definition 


The  state  for  the  cluster  was  initially  composed  of  the 
position  and  velocity  components  with  respect  to  the 
rotating  reference  point.  That  is,  until  Captain  Ward 
discovered  that  the  downrange,  y,  component  of  each 
satellite  was  unobservable.  The  component  was  removed  from 
satellite  #l's  state  and  for  satellites  #2 -s,  the  component 
was  replaced  with  a  relative  down-range  component  measured 
with  respect  to  satellite  #1.  The  state  vector  for 
satellite  #1  is  updated  to 


2  i 

7] 


(58) 


While  the  state  vectors  for  satellites  #2 -s  appear  as 


Ay, 

y< 


(59) 


Where  Ay,  is  determined  by  subtracting  the  y  solution  for 

satellites  #2-s  from  the  y  solution  for  satellite  #1. 
Therefore, 


Ay,(  +  )  =  y  ,(  +  )-y,(  +  ) 

=  [6(sinip-\ij)](xl -x,)  +  (yl(-)-y,(-))+  (60) 
"2 

-(cos  v -  1 ) 

Ln 


(*>  -*,)• 


4  3y 

-sin  - 

H  ti 


(y  i  -  y.) 


The  generalized  state  transition  matrix,  4>,  becomes 


Where 


= 


<t>  = 


+  i 

4> 


4> 


12 


13 


0 

<{>2 

0 


4>  J  s 


4  -  3cos  v 
0 

3r]sin  v 
6r|(cos  v  -  1 ) 
0 


0 

0 

4>. 


0  0 


0 

COS 

0 

0 

-  ti  si  n  \p 


si n  v 

TI 

0 

cosip 
■  2  sin  ip 
0 


0 

0 

0 


-( 1  -  COS  v) 

TI 

0 


2sin 

-  3+  4cosip 
0 


(61) 


0 

sin  ip 

h 

0 

0 

cos  y 


(62) 
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4-  3cos  ip 

0 

0 

sirup 

T] 

2 

-(  i  -  cos  ip) 

X] 

0 

-  6(sin  ip  -  ip) 

1 

0 

2 

—  (cosip  -  1  ) 

T] 

4  3  ip 

—  sin  ip  +  — 

n  n 

0 

0 

0 

cosip 

0 

0 

sin  ip 

n 

3r]sinip 

0 

0 

cos  ip 

2 sin  ip 

0 

6r|(cos  ip  -  1  ) 

0 

0 

-  2sin  ip 

-  3  +  4cos ip 

0 

0 

0 

-  r]  sin  ip 

0 

0 

cos  ip  _ 

(63) 

and 


0 

0 

0 

0 

o' 

6(sinip-ip) 

0 

2 

-(cosip  -  1  ) 

T] 

4 

-sin  ip 
r) 

3ip 

T| 

0 

+  u 

= 

0 

0 

0 

0 

0 

(64) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0_ 

Additionally,  the  h  vector  and  the  H  matrix  are  modified  to 


yield 


VC^i  -x2)2  +  Ay22  +  (2, -z2y 

= 

i 

.  .  ^r 

i 

i-xs)2  +  Aysz  +  (2,-2s)2_ 

i — 
Sr 

1 

1 _ 

H  ,  H  [  0  ...  0 
fi2  0  H2  ...  0 


H 0  0 


H 


s  —  1 


(65) 


(66) 
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where 


0  0  0 


(67) 


Xi-x^A  Ayit| 

h,  J  h, 


0  0  0 


(68) 


Once  minor  changes  are  made  to  the  truth  model  and  the 
estimator  algorithms,  the  computer  program  is  exactly  as 
it  was  when  Captain  Ward  generated  his  results . 
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III.  Observability  Analysis 


The  search  for  unobservable  states  was  complicated  by 
the  use  of  the  U-D  covariance  factorization  algorithm, 
since  the  filter  does  not  directly  invert  any  matrices. 
Therefore,  the  filter  was  replaced  by  a  non-linear  least 
squares  filter  to  take  advantage  of  the  fact  that  the 
matrix  ]TTtQ''T  (to  be  developed  shortly)  is  inverted 
during  the  estimation  process. 

3.1  Non-Linear  Least  Squares  Estimation 


The  equations  of  interest  in  the  non-linear  least 
squares  problem  are  the  update  to  the  reference  trajectory 
equation,  the  covariance  of  the  estimate  equation,  and  the 
residual  equation.  The  update  to  the  reference  trajectory 
is  given  by: 

6  x(t0)  =  (7TQ'17)'1TTQ'1r  (69) 

Where  for  convenience  the  matrix  product  H4>(t , ,  t0)  was 

redefined  as  7.  The  vector  r  is  the  residual  vector.  The 
trajectory  estimate  is  then  given  by: 

*(U  =  xr„,(U  +  6x(fJ  (70) 

The  estimate  covariance  is  given  by: 
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(71) 


=  C TTQ-'T ) 

Finally,  the  residual  vector  is  given  by  the  difference  of 
the  observed  and  the  calculated  data  vector. 

r  =  z  -  G(h(xa0). ',)•*.)  (72) 

The  variables  that  appear  on  the  right  hand  sides  of  the 
above  equations  were  previously  defined  in  Chapter  2.  The 
inversion  of  the  matrix  product  TtQ~xT  is  the  foundation  of 
the  calculations  that  appear  in  eqns  69-71.  If  the  matrix 
is  singular,  a  zero  eigenvalue  will  exist  (1:357),  and  the 
estimate  covariance  and  the  state  update  defined  by  eqns  69 
and  71,  respectively,  will  be  undefined.  In  other  words, 
unobservable  states  are  present  if  the  matrix  is  singular. 
Therefore,  the  ability  to  successfully  invert  the  matrix 
TtQ~xT  is  the  key  to  the  removal  of  the  unobservable 
states . 

3.2  Analysis 

The  analysis  initially  consisted  of  inverting  the 
matrix  ^TTQ"fT  using  a  Gaussian  elimination  with  maximal 
pivoting  algorithm.  Several  assumptions  were  made  to 
simplify  the  analysis:  namely,  the  constellation  consisted 
of  only  two  satellites  which  greatly  simplifies  the 
problem;  the  data  vector  consisted  not  of  the  range  between 
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the  satellites,  but  consisted  of  the  actual  components  of 
the  state  and  the  individual  data  measurements  were 
considered  independent  and  each  equally  contributing  to  the 
estimate  to  eliminate  any  questions  about  the  data  causing 
observability  problems;  and  finally,  the  initial  reference 
state  was  defined  as  the  initial  true  state  determined  by 
the  truth  model  thus  reducing  the  analysis  to  a  single 
iteration  of  the  non-linear  least  squares  estimation 
algorithm. 

The  first  attempt  to  evaluate  the  inverse  of  ]TTtQ~'T 

yielded  a  singular  matrix  as  expected,  confirming  earlier 
suspicions  of  state  unobservability.  The  next  step  was  to 
redefine  the  state  and  update  all  equations  defined  with 
respect  to  the  new  state  components  (this  task  follows  the 
steps  outlined  in  Section  2.3).  But  which  state  components 
should  be  removed?  An  examination  of  the  ]TTtQ'1T  matrix 
for  which  the  Gaussian  elimination  process  failed  yielded  a 
matrix  of  the  following  form. 
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0 

where  X  is  a  non-zero  value.  The  2  and  z  components  were 
removed  in  an  assumption  that  an  absolute  horizontal  plane 
of  reference  could  not  be  determined. 

The  new  state  consists  of  the  x,  x,  and  y  components 
for  satellites  #1  and  #2  and  the  relative  components 
Ay ,  Az, and  A z  between  satellites  #1  and  #2.  The  new 
equations  for  the  non-linear  least  squares  analysis  assume 
the  following  forms 


x 


h 


(74) 
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0 


(75) 


<£>  = 


^12  ^2 


4>.  = 


4  -  3cosip 


sin  ^ 


-( 1  -cosip) 

r\ 


T] 

3r|sinip  cosip  2sinip 
6r|(cosip-l)  -2sinip  -3+4cosipJ 


0 


0 


0 


2  4  3ip 

6(sinip-ip)  —  (cos  ip  —  1 )  -simp - 

T1  ti  T] 


0 

0 

0 

0 


0 

0 

0 

0 


0 

0 

0 

0 


(76) 


(77) 


4>, 


4-3cosip  0 

-  6(sin  ip  -  \p)  1 

0  0 


3 q sin  \p  0 
6r|(cos ip  -  1  )  0 

0 


0 

0 

cos  ip 

0 

0 


0  -  q  sin  ip 


simp 


—  (cos ip  -  1  ) 

n 

o 


cos  ip 
-  2sin  ip 
0 


-(  1  -  cos  ip) 

n 

4  3\p 

—  sin  ip  +  — 

Tl  q 

0 


2sin  ip 
-  3  +  4cos ip 
0 


0 

0 

sin  \p 

n 

0 

0 

cos  ip 
(78) 


Again,  the  inversion  of  ^TtQ~'T  failed,  yielding  a 


singular  matrix.  The  Gaussian  elimination  process  reduced 
the  matrix  to 
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(79) 
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000000000 
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000000000 
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One  may  examine  eqns  73  and  79  to  conclude  that  the  zero 

rows  corresponding  to  the  z2  and  z2  components  were 

removed  by  introducing  relative  2  and  2  components  into 

the  state.  Based  upon  this  conclusion,  the  zero  rows  for 

the  x2,X2,andy2  components  may  be  removed  by  introducing 

relative  A.v,  Ax.andAy  components  into  the  state.  This 

discovery  warranted  a  re-examination  of  the  solutions  to 

the  Clohessy-Wiltshire  equations.  If  the  solutions  for 

satellite  #1  and  #2  are  subtracted  from  each  other  and  if 

like  terms  are  collected,  one  obtains  the  following  set  of 

equations . 


Ax(  +  )  =  (4- 3cosip)A.v0(-)  +  — Ax„(-)  +  -(  1  -  cos\p) Ay 0(-)  (80) 

n  n 

Ay(  +  )  =  [6(sinv-v)]A.v0(-)  +  Ay0(-)  + 

2  f  4  3uA 

-(cosv-  1  )A.v0(-)  +  ^-sin  V  —  jAy0(-)  (81 ) 
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(82) 


ip 

Az(+)  =  costyAz0(-)  +  sin  -  Ai0(~) 

0 

Ax(+)  =  3r|Sin  tyAx0(-)  +  costy  Ax0(-)  +  2sin  tyAy0(-)  (83) 

Ay(  +  )  =  6t](cosij;  -  1  )Ax0(-)  -  2sin  Ax0(-)  +  (-3  +  4costy)Ay 0(-Q84) 

Az(  +  )  =  -rjsin  Az0(-)  +  cosipAz0(-)  (85) 

This  set  of  equations  indicates  that  there  is  no 
possibility  for  determining  the  position  and  velocity  state 
components  separately  for  satellites  #1  and  #2  based  upon  a 
range  measurement  between  the  two  satellites.  For  example, 
examining  the  Ax(+)  equation,  an  increase  in  the 
Xi.Xi.andy,  components  and  a  corresponding  increase  in  the 
x2,x2,andy2  components  would  yield  the  same  solution  as  if 
no  change  had  occurred  at  all.  In  other  words,  satellite 
#1,  the  filter's  host  satellite,  considers  itself  at  the 
origin  of  the  reference  frame  and  determines  the  relative 
motion  of  the  other  satellite  with  respect  to  itself.  Once 
identified,  the  discovery  seems  obvious  when  one  considers 
that  range  is  a  relative  measurement  and  by  itself  does  not 
yield  any  position  information. 

Equations  80-85  constitute  the  new  dynamics  model  that 
the  filter  will  use  during  the  estimation  process.  The  new 
state  vector  and  state  transition  matrix  for  a  cluster  of  s 
satellites  assume  the  following  forms 
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(86) 


(87) 
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(88) 
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Additionally,  the  h  vector  and  the  H  matrix  are  now  defined 


as 


(89) 


(90) 


(91) 
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IV.  Performance  Analysis 


Once  the  observable  states  are  defined  and  the  filter 
algorithm  completed,  the  next  step  is  to  tune  the  filter 
and  conduct  a  Monte  Carlo  simulation  in  order  to  assess  the 
filter's  performance.  The  tuning  of  the  filter  consists  of 
varying  the  diagonal  elements  of  the  dynamics  noise 
covariance  matrix  Qd‘,  a  rather  simplified  approach.  Once 
the  filter  is  adequately  tuned,  a  Monte  Carlo  simulation 
will  be  conducted.  The  constellation  size  will  be  varied 
and  so  will  the  random  number  seed  used  to  corrupt  the 
range  data.  Fifteen  test  cases  will  be  evaluated  for 
constellations  of  two,  five,  and  ten  satellites.  The 
positional  errors  between  the  truth  model  and  the  estimator 
will  be  determined  at  each  time  step.  Average  error  versus 
time  and  average  true  error  and  filter  covariance  versus 
time  plots  will  be  generated.  The  average  error  should  be 
near  zero  and  the  average  true  error  and  the  filter 
covariance  should  be  approximately  equal. 

The  filter's  performance  is  examined  by  comparing  the 
values  of  the  true  error  and  the  covariance  of  the 
estimate.  The  true  error  is  the  magnitude  of  the  position 
difference  between  the  truth  model  state  and  the  estimate 
state . 
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TRUE  ERROR  =  {e\  +  J\  ~  J\  (92) 

where  f,  is  the  difference  between  the  Ax,  Ay,  or  A 2 

component  of  the  true  state  and  the  estimate  state.  The 
estimate's  covariance  is  equal  to  the  square  root  of  the 
sum  of  the  squares  of  the  eigenvalues  associated  with  the 
position  components  of  the  estimated  state. 

0  =  \j  ( Eigenva  lue(x))2  +  ( EigenvalaQ^y ))2  +  (E  igenvalue(z))2 

(93) 


4 . 1  Filter  Tuning 

The  first  step  in  tuning  the  filter  is  to  assess  the 
filter's  performance  with  the  elements  of  the  dynamics 
noise  matrix  set  to  zero.  The  true  error  and  the 
covariance  are  plotted  together  in  Figure  2  for 
approximately  20  orbits  with  a  time  step  of  300  seconds. 
After  the  initial  transient,  the  covariance  exponentially 
decayed  approaching  zero.  The  estimate  approaches 
perfection  as  the  covariance  approaches  zero.  Once  the 
covariance  reaches  zero,  the  estimate  will  no  longer  change 
with  time  and  a  state  of  ignorance  will  exist  concerning 
the  future  behavior  of  the  estimated  state.  Therefore,  the 
filter  must  be  tuned  so  that  the  covariance  on  the  average 
will  neither  decay  to  zero  nor  grow  to  infinity  (7:84). 
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TIME  (SECONDS) 


Figure  2.  True  error  and  covariance  as  a  function 
of  time  with  zero  dynamics  noise. 


The  process  of  tuning  was  confined  to  the  diagonal 
elements  of  the  dynamics  noise  covariance  matrix  Qd  with  a 
constellation  consisting  of  two  satellites.  The  first  task 
was  to  establish  a  ballpark  initial  value  for  the  diagonal 
elements.  The  range  between  satellites  was  assumed  to  be 
accurate  to  within  one  centimeter.  Therefore,  one  may 
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determine  an  acceleration  value  that  will  yield 
approximately  a  one  centimeter  error  in  position  after  one 
orbit  by  solving  the  following  equation  for  a. 


0.01  meters 


1  ? 
-a(Period) 


(94) 


The  period  is  approximately  6300  seconds  and  the  calculated 
value  of  "a"  is  approximately  5X  10~10 A^M/SEC2.  Initially, 
only  the  diagonal  elements  corresponding  to  the  relative 
velocity  components  were  changed.  Therefore,  the  elements 
assumed  a  value  equal  to: 

Q  d  =  (a Afupda(e )2  (95) 


The  update  time  was  300  seconds  and  the  new  diagonal 
elements  were  approximately  2. 25X  10~I4AAf2/SEC2.  The 
exponential  decay  of  the  zero  dynamics  noise  filter  was 
removed,  but  the  true  error  curve  was  considerably  below 
the  covariance  curve.  Therefore,  the  filter  was 
overestimating  the  error  of  the  estimate  (see  Figure  3). 
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TRUE  ERROR 


TIME  (SECONDS) 


Figure  3.  True  error  and  covariance  as  a  function  of  time 
with  dynamics  noise  on  the  order  of  2.25A'  10‘14 A'AI2/SEC2 


Several  values  later,  the  filter  was  successfully  tuned 
with  an  acceleration  value  of  1  X  10'11  KM /SEC2  and  a  noise 
covariance  of  9X  10~18£A12/SEC2.  Figure  4  illustrates  the 
desired  result  that  the  true  error  and  the  covariance 
curves  overlap  one  another. 
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Figure  4 .  Filter  tuned  with  dynamics  noise  covariance  of 
9  A'  1 0" 18  A  Ai  2 /SEC2. 

The  diagonal  elements  of  the  noise  covariance  matrix 
corresponding  to  the  position  components  were  varied  but 
did  not  yield  any  appreciable  benefits.  For  all  subseguent 
testing,  these  elements  remained  set  to  zero. 

4.2  Monte  Carlo  Simulation 


Cnee  the  filter  was  successfully  tuned,  a  15-sample 
Monte  Carlo  simulation  was  conducted  for  constellations  of 
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2,  5,  and  10  satellites.  The  purpose  of  the  simulation  is 
to  test  the  filter's  performance  based  upon  varying  initial 
conditions.  The  truth  model  generated  the  true  state  and 
the  range  data  between  the  host  satellite  and  the  remaining 
satellites  at  each  time  step.  The  initial  estimator  state 
was  set  equal  to  the  initial  state  determined  from  the 
truth  model .  The  random  number  generator  seed  used  during 
the  corruption  of  the  range  data  was  changed  for  each 
simulation  run.  For  each  constellation  size  and  for  each 
simulation  run,  the  positional  errors  between  the  truth 
model  and  the  estimator,  and  the  filter  covariance  were 
calculated  and  stored  within  separate  data  files.  The 
final  results  consisted  of  an  average  error,  an  average 
true  error,  and  the  filter  covariance  at  each  time  step. 

For  constellations  of  5  and  10  satellites,  the  average 
errors  and  covariance  were  determined  for  three  different 
groupings  of  satellites.  For  example,  the  error  data  was 
calculated  between  satellites  #1  and  #2,  #1  and  #3,  and  #1 
and  #5  for  the  five  satellite  constellation. 

The  average  error  is  determined  by  summing  the 
positional  errors  from  the  15  data  files  for  each  time  step 
and  dividing  by  15. 

E  =  Z(£'x+  Ey +  f2),  (96) 
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The  average  true  error  is  determined  by  summing  the  square 
of  the  positional  errors  from  the  data  files,  dividing  by 
15  and  taking  the  square  root. 

The  covariance  calculated  by  the  estimator  at  each  time 
step  is  recorded  in  an  additional  data  file.  The  average 
error  should  be  approximately  equal  to  zero  and  the  average 
true  error  and  the  filter  covariance  should  be 
approximately  equal. 

4.2.1  Two  Satellite  Cluster 


Figures  5  and  6  depict  the  average  error  and  the 
average  true  error  and  covariance  results  for  a  two 
satellite  cluster.  As  expected  the  average  error  curve  is 
near  zero  and  the  average  true  error  and  the  filter 
covariance  curves  are  nearly  equal  as  a  function  of  time. 
The  average  true  error  approaches  a  value  of  approximately 
3  centimeters  or  three  times  the  range  measurement  error 
and  is  well  within  the  25  meter  accuracy  requirement 
(3:vii).  The  steady  state  filter  performance  exhibits  an 
oscillatory  nature  with  a  period  equal  to  the  orbital 
period  of  the  cluster  (approximately  6300  seconds). 
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Figure  5.  Average  error  versus  time  for  a  two  satellite 
cluster . 
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Figure  6.  Comparison  of  average  true  error  and  covariance 
versus  time  for  a  two  satellite  cluster. 


4.2.2  Five  Satellite  Cluster 


Figures  7-12  depict  the  average  error  and  the 
average  true  error  and  covariance  results  for  a  five 
satellite  cluster.  The  filter  was  not  re-tuned.  The 
diagonal  values  of  the  dynamics  noise  matrix  from  the 
'tuned'  two  satellite  constellation  case  were  assigned  to 
the  diagonal  elements  of  the  three  additional  satellites. 
Figures  7  and  8  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  two.  Figures  9 
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and  10  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  three.  Figures  11 
and  12  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  five.  The  filter 
performs  as  desired  with  an  average  true  error  on  the  order 
of  two  to  three  times  the  range  measurement  error.  The 
filter  continues  to  exhibit  the  oscillatory  behavior 
previously  noted. 
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Figure  7 .  Average  error  between  satellites  1  and  2 
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Figure  8.  Comparison  of  average  true  error  and  covariance 
versus  time  for  satellites  1  and  2. 
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Figure  9.  Average  error  between  satellites  1  and 


TEoooo 


3. 


45 


AVG  TRUE  ERROR 
COVARIANCE 


TIME  (SECONDS) 


Figure  10.  Comparison  of  average  true  error  and 
covariance  versus  time  for  satellites  1  and  3. 
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Figure  1.  Average  error  between  satellites  1  and 
5. 


Figure  12.  Comparison  of  average  true  error  and 
covariance  versus  time  for  satellites  1  and  5. 
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4.2.3  Ten  Satellite  Cluster 


Figures  13-18  depict  the  average  error  and  the 
average  true  error  and  covariance  results  for  a  ten 
satellite  cluster.  Again,  the  filter  was  not  re-tuned  and 
the  additional  diagonal  elements  of  the  dynamics  noise 
matrix  were  assigned  the  same  values  as  previously  used. 
Figures  13  and  14  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  two.  Figures  15 
and  16  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  five.  Figures  17 
and  18  are  the  error  plots  for  relative  position 
determination  between  satellites  one  and  ten.  The  expected 
data  trends  are  obtained  for  the  ten  satellite  cluster 
case.  The  average  true  errors  are  approximately  two  to 
three  times  the  range  measurement  error  and  again  the 
oscillatory  behavior  of  the  filter  is  exhibited  as 
previously  discussed. 
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Figure  14.  Comparison  of  average  true  error  and 
covariance  versus  time  for  satellites  1  and  2. 
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Figure  15.  Average  error  between  satellites  1  and 
5. 
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Comparison  of  average  true  error  and 
versus  time  for  satellites  1  and  5. 
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Figure  17.  Average  error  between  satellites  1  and 

10. 
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Figure  18.  Comparison  of  average  true  error  and 
covariance  versus  time  for  satellites  1  and  10. 
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V.  Conclusions  and  Recommendations 


Relative  position  determination  of  a  satellite 
cluster  using  an  on-board  estimator  is  possible  and  yields 
good  results  tested  against  a  two-body  astrodynamics  model. 
The  state  of  the  cluster  can  only  involve  relative  position 
and  velocity  components  between  the  filter  host  satellite 
and  the  remaining  satellites.  The  U-D  Covariance 
Factorization  Kalman  Filter  was  tuned  and  subjected  to  a 
15-sample  Monte  Carlo  Simulation.  The  error  results 
presented  illustrate  good  behavior  of  the  average  error, 
the  average  true  error,  and  the  filter  covariance.  The 
average  true  errors  were  approximately  three  times  the 
range  measurement  error. 

Further  investigations  into  this  topic  should 
include  a  full  tuning  analysis,  a  filter  robustness  test,  a 
test  of  the  filter's  performance  using  a  more  accurate 
truth  model  and  an  investigation  into  the  filter's  cyclic 
behavior.  The  variation  of  non-diagonal  components  of  the 
dynamics  noise  covariance  matrix  should  be  investigated  for 
any  significant  gain  in  performance  or  accuracy.  A  full 
test  of  the  filter's  robustness  should  be  conducted  to 
analyze  the  filter's  performance  with  a  poor  initial  guess 
for  the  initial  state,  and  with  periods  of  highly 
inaccurate  range  data.  A  test  with  an  extremely  accurate 
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perturbations  model  for  the  truth  model  should  be  conducted 
as  a  preliminary  test  to  establish  the  filter's 
flight-readiness.  Finally,  the  filter's  cyclic  behavior 
should  be  investigated  to  fully  understand  the  filter's 
performance . 
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